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ANALYSIS  OF  THE  BINARY  EUCLIDEAN  ALGORITHM 


Richard  P.  Brent 
Australian  National  University 
and  Carnegie -Mellon  University 


1 . Introduction 

The  binary  Euclidean  algorithm  of  Silver  and  Terzian 
[62]  and  Stein  [67]  finds  the  greatest  common  divisor  (GCD) 
of  two  integers,  using  the  arithmetic  operations  of  subtrac- 
tion and  right  shifting  (i.e.,  division  by  2).  Unlike  the 
classical  Euclidean  algorithm,  nc  divisions  are  required. 
Thus,  an  iteration  of  the  binary  algorithm  is  faster  than  an 
Iteration  of  the  classical  algorithm  on  many  binary  computers. 

The  classical  algorithm  has  been  exhaustively  analyzed 
from  the  time  of  Gauss:  see,  for  example,  Dixon  [70,  71], 

Gauss  [12],  Heilbronn  [68],  Khinchin  [35a,  35b,  36],  Kusmln 
[28],  Levy  [29],  Szllsz  [61],  Tonkov  [74]  and  Wirsing  [74], 

A good  survey  is  given  in  Knuth  [69],  The  theory  of  the 
binary  algorithm  is  much  ess  satisfactory.  Knuth  [69]  ana- 
lyzed a "lattice-point"  model  which  is,  unfortunately,  only 
a crude  and  pessimistic  approximation  to  the  actual  algorithm. 
In  this  paper  we  analyze  a continuous  model  of  the  binary  al- 
gorithm and  find  the  expected  number  of  iterations.  The  re- 
sults agree  with  the  observed  behavior  of  the  algorithm  much 
better  than  those  predicted  by  Knuth' s "lattice-point"  model. 

The  binary  Euclidean  algorithm  for  finding  the  GCD  of 
positive  integers  u and  v is  given  in  Knuth  [69,  Sec.  4.5.2, 
Alg.  B],  After  steps  Bl  to  B5  of  the  algorithm  have  been 
performed  once,  the  problem  is  reduced  to  that  of  finding 
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the  GCD  of  two  odd  Integers.  Thus,  we  assume  here  that  u 
and  v are  odd,  and  the  algorithm  is  as  follows. 

RS  Binary  Algorithm 
[n  - 0; ] 

LI : t *-  |u  - v| ; 

if  t ■ 0 then  return  u as  the  GCD  and  halt; 

L2 : t *-  t/2; 

if  t is  even  then  go  to  L2; 

L3:  [n  ♦"  n + 1 ; ] 

if  u 2 v then  u «-  t else  v *-  t; 
go  to  LI . 

The  statements  in  square  brackets  are  not  essential.  We  say 
that  one  "iteration"  is  one  execution  of  step  L3,  so  n counts 
the  number  of  iterations.  To  distinguish  the  different  val- 
ues taken  by  the  variables  u and  v,  we  let  uq  be  the  value  of 
u at  Iteration  n,  etc.  Step  L2  is  executed  twice  as  often  as 
step  L3,  on  the  average,  but  the  L2  loop  merely  shifts  t 
right  until  it  becomes  odd,  and  this  may  be  done  efficiently 
on  a binary  computer. 

Let  x ■ min(u  , v )/max(u  , v ),  and  let  F (x)  be  the 
n n n n n n 

probability  distribution  function  of  x^.  We  assume  that  u^ 
and  Vq  are  uniformly  and  independently  distributed  in  (0,  N) 
(with  the  constraint  that  they  are  odd),  and  consider  the 
continuous  approximation  obtained  by  letting  N -*  “.  In  Sec- 
tion 2 we  derive  a recurrence  relation  for  the  continuous 

distributions  F (x). 

n ^ 

In  Section  3 we  show  that  F (x)  » OL(x)lg(x)  + 0 (x), 

n n n 

where  a (x)  and  0 (x)  are  analytic  and  satisfy  certain  recur- 
n n 

rence  relations.  An  explicit  expression  for  a (x)  is  given 

n 

* 

Throughout  this  paper,  lg(x)  denotes  iog2(x). 


PUf,  I.  I,  im.reyiwywrry* 


r"  !'^v,f' i,M  ,,'T^r?'  * 
f 


In  Section  4. 

In  Section  5 we  consider  the  equivalent  recurrence 
f^l  ■ Tf^  for  the  probability  density  functions 

f^(x)  ■ F^(x).  We  8^ow  ll£n«-£„+lH<  l^+r£nll£or* 

certain  norm.  Numerical  evidence , described  in  Section  7, 

suggests  that  convergence  is  rapid.  Thus,  it  is  likely  that 

f tends  to  a limiting  density  f , though  we  have  not  been 
n 50 

able  to  prove  this. 

The  expected  number  of  iterations  is  asymptotically 
Klg(N)  for  large  N,  and  an  expression  for  the  constant  K is 
derived  in  Section  6.  The  theoretical  value  of  K ~ 0.706 
agrees  with  values  obtained  numerically  for  moderate  values 
of  N.  The  numerical  results  are  described  in  Section  7. 

Finally,  in  Section  8 we  consider  another  algorithm 
which  uses  only  shifts  and  subtractions.  The  algorithm  uses 
left  shifts  (i.e.,  multiplication  by  2)  instead  of  right 
shifts,  so  we  call  it  the  left-shift  binary  Euclidean  algo- 
rithm (LS  algorithm  for  short).  We  show  that  the  expected 
number  of  Iterations  is  slightly  greater  than  for  the  (right- 
shift)  binary  Euclidean  algorithm.  However,  the  LS  algorithm 
is  worth  considering  for  use  on  a computer  with  a "normalize" 
Instruction,  as  the  left-shifting  loop  may  be  replaced  by  one 
Instruction.  Either  of  the  binary  algorithms  could  be  imple- 
mented in  hardware  (or  microprogrammed)  with  approximately 
the  same  expense  as  integer  division. 

We  consider  only  single-precision  integer  GCD  computa- 
tions here.  For  polynomial  and  multiple-precision  Integer 
GCD  algorithms,  see  Collins  [74],  Schonhage  [71]  and  Knuth 
[69]. 

2.  The  Recurrence  for  F 
n 

For  notational  simplicity  we  write  u for  ufl  and  u’  for 


— 
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u . , , etc.  Also,  there  Is  no  loss  of  generality  in  assuming 
n+1 

that  u s v.  The  iteration  terminates  if  u - v,  so  we  assume 

that  u > v.  Thus,  x « v/u,  t ■ 2 k(u~v),  and 

x'  ■ min(t,v)/max(t,v) , where  k 2 1 is  chosen  so  that  t is  an 

odd  integer. 

Let  P(E)  denote  the  probability  of  an  event  E.  By  defi- 
nition, (y)  - P(x’  £ y),  but  x'  - min(t/v,  v/t),  so 

- P(t/v  s y V v/t  £ y) 

■ P(t  ivy  V v £ ty) . 
shown  that,  for  K - 1,2,..., 

- K)  - 2‘K. 


■ ^ 2 k P(2  k(u-v)  £ vy  V v £ 2 k(u-v)y) 
k-1 

00 

■ ^ 2”kP(2’k(1 -x)  £ xy  V x £ 2“k(l-x)y). 
k-1 

Since  x € (0,1),  we  have  2 k(1-x)  £ xy  iff  x 2 l/(l+2ky),  and 
x £ 2 k(1-x)y  iff  x £ l/(l+2k/y),  Also,  assuming  y £ (0,1), 
we  have  l/(l+2k/y)  £ l/(l+2ky).  Thus,  from  (2.5), 

00 

(2.6)  F^y)  - ^2'k[l-P(l/(l+2k/y)  £ x £ l/(l+2ky))]. 
k-1 

Since  x has  distribution  function  F^,  this  gives  the  interest- 
ing recurrence  relation 


irr  i 


(2.1)  F^W 


(2.2) 

It  may  be 

(2.3)  lim  P(k 
N—> 

Thus, 

(2.4)  Fttfl(y) 


(2.5) 


v'"'1  «pi  yri,.  » 
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^‘Q'"  '• 

for  n 3s  0 and  y 6 [0,1  ]. 

The  corresponding  recurrence  for  the  classical  algorithm 
la 

00 

(2.8)  C^Cx)  - £ (G„OA)  - Gn(l/(k+x»]. 

k-1 

This  was  derived  by  Gauss  [12],  who  conjectured  that 

(2.9)  lira  Gn(x)  - lgd+x)» 

n*4® 

ohlch  <u  proved  by  Kusmln  [28).  Sharper  results  »ere  later 
obtained  by  u'vy  [29]  and  Stilsx  [61].  Finally.  Ulrsing  [74] 
proved  that 

(2.10)  Gn(x)  - lg(1+x)  + 0(Xnx(1-x)) 

as  n uniformly  for  all  x € [0,1],  where  X ~ 0.3036630029 

is  a certain  constant  in  (0,1). 

We  conjecture  that  a similar  result  holds  for  Fr(x) . 

For  a reason  which  will  be  clear  later,  the  term  x(l-x)  in 
(2.10)  must  be  replaced  by  x|ln(x)|. 

Conjecture  2.1 

There  exists  Fb(x)  ■ lim  Ffl(x) , and 

n-** 

(2.11)  Ffl(x)  ■ F^x)  + 0(Xnx  | ln(x)  | ) 

uniformly  for  all  x € (0,1],  where  X is  some 


(2.7) 


r„+i(y)  ' 1 


r ■ v 
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constant  In  (0,1).  . 

Th.  theoretical  evidence  tor  Conjecture  2.1  Is  given 

the  next  three  section.,  end  some  numerical  evidence  1.  given 
in  Section  7 . 

Differentiating  (2.7),  «.  obtain  th.  recurrence 


(2.12)) 


£n+l(X>  ’ L 
k-1 

f0(«)  - 1 


fc)2  £46 + (i)  c;o] 


for  the  probability  density  functions  fn(x)  - F’(x), 

X € (0  1],  n * o.  The  recurrences  (2.7)  and  (2.12)  are 
equivalent,  but  In  Section  3 »e  prefer  to  work  with  (2.7)  and 
consider  the  four,  of  Fn(x).  Results  for  £„<*>  “slly  e‘ 


duced  by  differentiation. 

3 1 The  Distribution  Functions  F 


The  following  theorem  gives  the  form  of  Fn(x)  for  finite 


Theorem  3.1 

For  all  n 2;  0 and  x £ (0,1  ]» 

(3.1)  Fn(x)  - <*n(x)lg(x)  + 8n(x>* 

where  on(x)  and  0n(x)  ere  analytic  and  regular  In  |x|  < 1. 

and  o (0)  - 0(0)  ' »•  Al5°-  “o(x>  ' ° anli 

n n 

(3.2)  2V](2x)  - o^M  - o„(^)  * 3fn(1)x- 
Proof 

Define  DQ(x)  - 0 and 


'W  * 
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(3.3)  D^Cx)  - I*'kr„Br)  • 

k»l  V X 

We  assume  that 

(3.4)  F (x)  ■ a (x)lg(x)  + B (x) , 

mm  m 

(3.5)  D (x)  - 1 + v (jcHg(>c)  + 6 (x) , 

m m id 

(3.6)  D (l/x)  - e (x)lg(x)  + 0 (x) , 

id  m m 

and 

<3‘7>  Fm(i^)-1+Vx) 

for  m < n,  where  om(x) , . . . ^(x)  are  analytic  and  regular  for 
|x|  < 1,  and  vanish  at  x ■ 0,  We  shall  prove  the  correspond- 
ing result  for  m ■ n,  so  (3.1)  will  follow  by  induction.  The 
results  (3.4)  to  (3.7)  are  trivially  true  for  m ■ 0,  so  we 
may  assume  n > 0. 

From  (2.7)  and  (3,3)  we  have 

(3.8)  F (x)  - 1 + D (l/x)  - D (x) , 

n n n 

jo  if  a (x) 0 (x)  are  regular  at  x - 0 we  must  have 

n n 

(3.9)  «n(x)  - eR(x)  - Yn(x) 
and 

(3.10)  B (x)  - d <x)  - 6 (x) . 

n n n 

From  (3,3)  we  also  have 

(3-">  2Bn(£)  - DJj)’  *■.-.(£)  • 

so  in  the  same  way  we  find  that 
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(3.12)  2<n(2x)  - ' V.  (ife) 

and 

(3.13)  2c  (2x)  + 2^  (2x)  - *5  (x) 

n n n 

By  the  Inductive  hypothesis,  the  right  si'.es  of  (3.12)  and 

(3.13)  are  analytic  and  regular  at  x ■ 0.  Let  the  Taylor 
expansion  of  a (x)  be 


(3. 14)  or  (x) 
m 


L 

j-i 


f 


and  similarly  for  ^^(x) ,. . . ,T^(x) . By  equating  coefficients 
we  see  that  analytic  solutions  ^(x)  and  ^n(x)  satisfying 
(3.12)  and  (3.13)  exLst,  and  are  given  by 


(3. 15)  e 

OtJ 


and 

(3 • 1 6)  6 


where  J ■ 1,2,...  . Thus,  en(*)  *nd  $n(x)  are  determined  by 

a . (x)  and  0 . (x) , and  are  analytic  and  regular  in  Ixl  < 1. 

n- i 0“ i 

From  (3.3)  and  (3.8), 
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(3.17)  Fn(y>  - 1 - l Fn_,(^)  + 1 

- 1 v*> ♦ f »J*)  • 

Substituting  y - l/(1+x)  gives 

<3-'8>  F„(t^)  - ' - 2 VlIjS)  + F Fn-1  (sk) 

- F DXife)  + F Dn(2+2x)- 

By  the  inductive  hypothesis, 

<3J9)  Fn-(H  * “n-lCH1^)  + 

f 2x\  f -2x  \ 

so  substituting  y "1 3(3+xjJ  andU(3+2x)J  8tveS  P°Wer  8erl®8 
for  Fn-l(j£)  ind  F»-1  respectively.  Also, 

<3-2°>  Dn(r+jr)  ■ *PtKt) + <¥) 

end 

(3.21)  Dn(2+2»)  - -«B(^  l»«*x)  ■ 

Ih“'-  F/tB  ' 1 + VK,> 


where  7^(x)  is  analytic  and  regular  in  |x|  <1. 

It  remains  to  consider  y (x)  and  6 (x).  From  (3.3), 

n n 

(5.22)  2o(f)-  Dn(x)  - 1 +!!„.,  (x), 

80 

(3.23)  2V[(j)  ' Vn(*>  - 0 


ill'V*UIRHKI 
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and 


(3.24)  26/f)  - 6n(x)  -2V/f)- V,(x), 
Thus,  we  have  the  analytic  solutions 


(3.25)  V.M  - V,.,*- 

and 

1 


(3.20)  6 


n.J 


Vi.j 


for  12  2.  The  constant  6 , may  be  determined  from  the  re* 

n»* 

la  t ions  0 . - 0 , - 5 , and 

ii|  I n,  i n y l 


(3.27)  ?l2)  " 1 - J Fn-lC2 ) “ 4 Fn-lCs  ) + 4 Dn(2)’ 


obtained  from  (3.10)  and  (3.17)  respectively. 

We  have  now  proved  (3.4)  to  (3.7)  for  m • n,  so  the 
first  part  of  the  theorem  follows  by  induction.  (3.2)  fol- 
lows easily  from  (3.9),  (3.12)  and  (3.23),  so  the  proof  is 
complete. 

It  is  interesting  to  obtain  an  explicit  formula  for 


F^(x).  First  we  need  a lemma. 


i&IBP*  LI 

If 


(3.28) 

D) 

(x) 

" L2 

-k/(H2kx) 

1 

k"1 

then 

(3.29) 

Di 

(X) 

■ xlgx 

+ ) + | - 

2 

x 

1+x  ' 

r 

L 

J-2 

for  0 < 

X 

< 2, 

, and 

4.  • ~ 


•WWipililppw 


W rp-fjv-i'T-TJI  1 1 WF^pp^T  ,r: 


H*r«l*r  <■■»•  «-r-n  «<>.*.* ,- 
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(3.30)  D^l/x)  - - £ *j$— 

j-1  2 -1 

for  |x|  < 2. 

Proof 

From  (3.5)  and  (3.6),  we  have 

(3.31)  (x)  - 1 + v^xHgtx)  + 5,(x) 

and 

(3.32)  D^l/x)  - «1(x)lg(x)  + (x) . 

Since  o0(x)  - 0 and  0Q(x)  - x,  (3.15)  gives  c^x)  - 0,  and 
(3.16)  gives  - (-1  )^f1/(2^+1 -1) • This  establishes  (3.30) 
From  (3.25),  •y]  (x)  - x.  Also,  since  T)Q(x)  - l/(1+x), 
(3.26)  gives 

(3.33)  6,  , “ (-1)J/(21“j-l) 

1 » J / 

for  J 2 2.  Thus 

r (-x)J 

(3.34)  D1  (x)  - xlgx  + 1 + 6,  ,x  - £ . 

J-2  ''2 

The  series  in  £3.34)  converges  for  |x|  < 1.  Subtracting  and 

adding  - ),  (-x)J  gives 

j-2 

x2  V ( ~x)  ^ 

(3.35)  Dt(x)  - xlgx  + 1 + 6ulx  - ^ “ 2.  T^l”  * 

j-2  “ 

where  the  last  series  converges  for  |x|  < 2.  By  analytic 
continuation,  (3.35)  holds  for  0 < x < 2.  The  constant 


f 

1 


HUM 


A. 
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6 • ^ may  be  determined  by  equating  (3.30)  and  (3.35)  with 

1 # * * 

x ■ 1.  Thus,  (3.29)  follows. 


Corollary  3.2 


F1  (x) 


-xlg(x)  + 


x(5x-1 ) 
6(1+x) 


+ 3 


I 


(-x)j  2j_1 


(2j“1-1)(2J+1-1) 


Proof 

This  follows  from  (3.8)  and  Lemma  3.1. 

In  principle  we  could  obtain  F^(x),  F^(x),  etc.  in  the 
same  way  as  F^ (x) . However,  the  details  become  very  compli- 
cated. The  situation  is  similar  for  the  classical  algorithm: 
see  Knuth  [69], 

Corollary  3.3 

For  all  n 2 0 and  some  x € [0,1],  F^j  t Fn(x)» 

Proof 

Suppose,  by  way  of  contradiction,  that  F^  (x)  ■ FR(x) 

for  all  x £ [0,1],  From  Corollary  3.2,  n f 0.  From  Theorem 

3.1,  or  (x)  - a (x).  Thus,  from  (3.2), 
ttr  I n 

<3-36>  ■ “n(i^  • 3fn(')11 

for  | x | < 1. 

Substituting  y ■ x/(l+x)  we  obtain 

(3.37)  an(y)  - 3n-1(y)  - 3(fn(D  - fn_, <i»y/0 -y) 

for  | y J < j.  By  analytic  continuation,  (3.37)  holds  for 

) y J < 1.  However,  from  (3.2)  it  follows  that  <*n(y)  *nd 

a , (y)  art  regular  at  y ■ 1,  so  we  must  have  f (1)  ■ f ,(1), 
n—  i n n— i 


fin 


'WWf 
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and  thua  a (y)  ■ a Ay)  • Continuing  in  this  way,  we  finally 
n n-1 


obtain  o^(x)  ■ aQ(x),  which  contradicts  Corollary  3,2 


(a  (x)  ■ x,  Oq(x)  ■ 0).  Thus,  the  original  assumption  was 


false, and  (x)  i Fn(x)  for  some  x € [0,1]. 


4.  Solution  of  the  Recurrence  for  a 


In  this  section  we  solve  the  recurrence  (3.2)  explicitly. 
The  method  used  here  can  obviously  be  generalized.  However, 


we  have  not  been  able  to  solve  the  recurrence  for  @ (x) 

n 


analytically. 

Define  p(0)  m 0 • 

p(2n)  - p(n)t 

and  p(2nfl)  ■ p(n)  J-  1. 

Thus,  p(n)  is  the  number  of  one-bits  in  the  binary  representa- 
tion of  n 2 0. 


Theorem  4.1 

Suppose  £>q(x)  • 0 and 


(4.1)  2cn+,(2x)  - V,(x)  - ajjg)  + cn+1x 


for  n 2 0,  where  c^Cg,...  are  constants,  cQ  - c_j  ■ 
and  a . , (x)  is  analytic  and  regular  at  x ■ 0.  Then 

Ctrl 


- o, 


2k-l 


(4.2)  o(x)-5.  J 2-k;[  -pui 

k-0  j-0  2+J* 


for  all  n 2 0 and  all  x 4 (-*,  -1], 


Note 


(4.1)  is  the  same  as  (3.2)  if  c . , ■ -3f  (1)  for  n 2 0. 


‘n+1  n ' 

Thus,  (4.2)  gives  an  explicit  solution  of  (3.2)  in  terms  of 


f0(1),  f1(1),...,fn_,0) 


•’  -riiiifiawMIgSMMiirtatMfc. 


n.'T" 


14 


Proof  of  Theorem  4.1 

The  result  ic  true  for  n - 0,  and  the  analytic  solution 
of  (4.1)  which  is  regular  at  x ■ 0 is  clearly  unique.  Thus, 
it  is  sufficient  to  verify  that  if  a (x)  and  (x)  are  de- 
fined by  (4.2)  then  (4.1)  holds.  From  (4.2)  we  have 


2Vi (2x)  ' Vm  m - 


2k+1_i 


x -k  r 

4 L 2 L 

k--1  j-0 


’n-H-p(J) 

2k+jx 


• 2k-1 

x \ -k  \ 

“ 4 L 2 c\ 

k-0  J-0 


2k+jx 


2k+1-1 


x \ -k  \ 

4 L. 


k-0 


j-2 


/ Cn4-1  -P(J) 

.tic  -k. 


2 +Jx 


Vix’ 


since  p(2  +J) 
follows. 


y 

p(j)  + 1 for  0 S J < 2‘,  Thus,  the  result 


g9r9Harg  4J 

Suppose  lim  f (1)  - f^O)  exists. 
n-*«  n 

iio  or  (x)  - a (x)  exists,  and 
n ® 

n-*«® 


Then 


-31(1) 


(4.3)  cy^x)  - y 


♦(x). 


where 
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For  simplicity  we  assume  x is  reel  end  positive,  though 
e similer  proof  goes  through  for  complex  x 4 (-»,  -1].  From 

(4.9)  we  heve 

OB 

(4.10)  |V,(x>  l 2'\.k 


Given  c > 0,  there  exists  m such  that  d s e.  Thus, 

m 

for  n i max(m,  m+lg(dQ/e)),  we  heve 

od  n-m  • 

r -it  r -v  \ _w 

/2Kdtfi/2Kd.+  / 2 o . 

n-k  ^ n-k  L,  n-k 

k"0  k*0  k*»n-m+1 

< 2c  + 2m“nd0  s 3c 

Thus,  lim  a (x)  exists,  and  the  limit  is  given  by  (4.3)  and 
n 

(4.4). n 

The  recurrence  (4.5)  may  be  verified  as  in  the  proof  of 
Theorem  4.1,  and  equating  coefficients  gives  (4.6).  Also, 
substituting 


(4.11)  ”7” — - 2*k  ^ (-2’kjx)n 

2 +JX 


in  (4.4)  and  equating  coefficients  gives  (for  n > 1) 
• 2*-l 

(4.12)  l 2**"  l Jn'\ 

k-1  j-1 

so  (4.7)  follows  from  ex.  1.2.11.2.4  of  Knuth  [68]. 
Corollary  4.2 

Suppose  lim  f (1)  ■ f (1)  exists,  and  that 

..  — n • 


w 
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(4.13)  fn(1)  ■ f.d)  + 0(Xn) 

as  n -*  ®,  where  X € (j»  1).  Then 

(4.14)  a (x)  ■ o (x)  0(Xnx) 

n ® 

and 

(4.15)  o'  (x)  - «'  (x)  + 0(Xn) 

n ® 

as  n -*  ®,  uniformly  for  all  x g [0,1]. 

Proof 

From  (4.10), 

eo 

I“nf1(x)  “ a«(x)l  “ 0(X“x)  l(2X)’k» 

k*0 

and  2X  > 1 , so  the  last  series  is  convergent.  The  proof  of 

(4.15)  is  similar. 

5.  Some  Convergence  Results 

We  define  a linear  operator  T,  mapping  the  Banach  space 
Lj(0,1)  into  itself,  by 

CO 

(5.1)  Tf (x)  - £ 
k-1 

Thus,  (2.12)  is 

<5*2)  fr*l-Tfn‘ 

We  write  f 2 0 if  f(x)  2 0 for  almost  all  x € [0,1]  (in 
the  sense  of  Lebesgue  measure).  Note  that  T is  a positive 
operator,  i.e.,  Tf  2 0 whenever  f 2 0. 


iUithtttf* 


This  proves  (5.3).  To  prove  (5.4) , we  merely  note  thst  ell 
the  inequalities  In  the  proof  of  (5.3)  become  equalities  if 


f 2 0. 


Proof 


This  is  immediate  from  Theorem  5.1  and  the  definition  of 


We  would  like  to  prove  that  the  iteration  (5.2)  con- 
verges to  a fixed-point  of  T.  Unfortunately,  the  theorems  of 
Schauder  (see  Simmons  [63])  and  Krein  and  Rutman  [43]  are  not 
applicable,  because  [f  £ (0,1)  | ||f  ||  ■ 1}  is  not  compact. 

Thus,  we  have  only  been  able  to  prove  the  weaker  result  given 
in  Corollary  5.2. 


Theorem  5.2 


Suppose  rhat  f is  continuous  on  (0,1),  changes  sign  at 
least  once,  does  not  vanish  on  any  finite  subinterval  of 
(0,1),  and  there  exists  e > 0 such  that  f(x)  ■ 0 has  no  solu- 
tion x £ (0,e],  Then 


(5.6)  Ifrf||  < ||f ||. 


Proof 


Suppose,  by  way  of  contradiction,  that  ||lf||  M||f jj.  Thus, 
all  Inequalities  in  the  proof  of  Theorem  5.1  must  be  equali- 
ties. Hence,  for  all  k 2 1 and  all  x £ (0,1),  we  have 


(5.7) 


<SX^) 


By  assumption,  f(x)  changes  sign  at  some  point  <p  £ (0,1) 
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There  exists  K 2 1 such  thet  cp  > . Suppose  k 2 K,  so 

1+2* 


qj  > 7~  . Then  there  exists  x.  £ (0,1)  satisfying 

1+2K  K 

. Thus,  from  (5.7),  f must  also  change  sign 
* k 

V < 2 . Since  k may  be  arbitrarily  large, 

this  contradicts  the  hypotheses  of  the  theorem.  Thus,  (5.6) 
must  hold. 


<p  - l/(l+2kxk) 

« ■ V<2“ 


figJPllgrYJi?. 

Let  e ■ f , , - f . Then 
n n+1  n 


(5.8)  IK.rt-,  II  < IM 


for  all  n 2 0. 


Proof 

From  (5.2),  e , , * Te  , so  we  have  only  to  show  that  e 
n+ 1 n n 

satisfies  the  conditions  of  Theorem  5.2,  From  Theorem  3.1, 

e (x)  " 6 (x)lg(x)  + £ (x),  where  a (x)  and  i (x)  are  analyt- 
n n n n n 

ic.  Also,  from  Corollary  3.3,  e^(x)  does  not  vanish  iden- 
tically. Thus,  e^(x)  is  continuous  on  (0,1)  and  does  not 
vanish  on  any  finite  subinterval  of  (0,1). 

Since 


(5.9)  J1  6 (x)dx  - J1  f (x)dx  - J°  f (x)dx  - 0 
0 n 0 n 0 n 

but  lie  |j  > 0,  e (x)  must  change  sign  at  least  once  on  (0,1). 
11  nM  n 

Finally,  from  Theorem  3.1  we  see  that  eR(x)  has  constant 
sign  on  (0,e],  for  some  sufficiently  small  e > 0.  Thus,  the 
conditions  of  Theorem  5.2  are  satisfied,  and  the  result  fol- 
lows. 

From  numerical  evidence  we  conjecture  that 


w.iii 
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\ 

* 


/ 


v 


j 

I 

I 

I 


f 


I 

I' 

j 

i 


(5.10)  IK*,  II  * X||»nll 

for  some  X € (0,1).  Unfortunately,  Corollary  5.2  does  not 
imply  (5.10).  If  (5.10)  is  true  then  (f^)  is  a Cauchy  se- 
quence and  the  limit  f exists. 

CO 

Corollary  5.3 

For  all  n i 20,  and  all  x £ [0,1], 

(5.11)  |F,*,(x)  - Fn(x)|  S |KJ|<  !0‘10. 

Proof 

|F,i(x)  - F (x)  | - IJ*  e (y)dy|s  |je  ||,  but  numerical  re- 
irT  i n WQ  n n *10 

suits  (described  in  Section  7)  show  that  I^qII  <10  , so 

the  result  follows  from  Corollary  5.2. 

From  now  on  we  assume  that  the  limiting  distribution 
F(b(x)  exists.  In  view  of  Corollary  5.3,  we  may  use  F2q(x) 
Instead  of  F (x)  for  all  practical  purposes. 

GO 


6.  The  Expected  Number  of  Iterations 

We  use  the  notation  of  Section  2.  Let  s ■ u+v  and 
a'  ■ u'+v' . Note  that 


(6.1)  a/s' 


(u+v)/ (u’+v1) 


Since  k 2 1,  s/s1  2 2,  so  the  maximum  number  of  iterations  is 
at  most  Ug(N)J.  The  example  u ■ 2m-1,  v ■ 1 shows  that  this 
bound  is  attainable.  For  another  example  see  Knuth  [69], 
exs.  4.5.2.27-28. 

Let  En  be  the  expected  value  of  ln(s/s').  From  (6.1), 


— 1 

I 1 
l ' 


! 

j 

I 

\ 

i 


f 


mi  inr 


J 
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so 

(6.2) 

E - ln2  f 
n 

where 

00 

(6.3) 

$(x)  - ^ 
kP2 

2 (1+x)  * 


n-1 

The  expected  value  of  ln(Sg/sn)  Is  ^ Ej.  Thus,  assum- 

j-0 

ing  the  existence  of  E ■ lim  E,,  the  expected  number  of 

j-®  J 

Iterations  for  odd  integers  u^,  S N is  asymptotically 
K lg(n)  is  N •*  «,  where 

(6.4)  K - ln(2)/E  . 


Approximating  E^  by  E^Q  and  evaluating  the  integral  in 
(6.2)  numerically  gives 

(6.5)  K ~ 0.705971246102. 


In  the  next  section  we  give  some  numerical  evidence  which 

suggests  that  the  expected  number  of  iterations  is 

K lg(n)  + 0(1).  This  is  not  surprising  if  Conjecture  2.1 

holds,  for  then  E ■ E + 0(X  n) . 
n “ 

7.  Numerical  Results 

The  recurrence  relation  (2.7)  was  solved  numerically  by 
three  different  methods.  All  computations  were  performed  on 
a Univac  1108  using  double-precision  (60-bit  fraction),  and 
the  numerical  results  given  by  the  different  methods  agreed 
to  the  accuracy  expected. 


k,.. 


. ; it'.  ‘M.- XVrifailitSiafaa*  • . i 


; 


■j 

} 


k 


1 

3 

'i 


| 
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A.  The  Recursive  Method 

This  is  the  most  obvious  method.  F (x)  is  evaluated  re- 

n 

cursively,  using  the  recurrence  (2.7)  with  the  Infinite  sums 
truncated  after  the  terms  become  negligible.  The  method  is 
only  useful  for  small  n,  as  the  computation  time  increases 
exponentially  with  n. 

B.  The  Discretization  Method 

If  F (x)  is  known  ".t  a finite  set  of  points,  say 
n 

x.  ■ 0 < x,  < x-  < ...  < x "I,  then  we  can  use  the  recur- 
0 1 Z m 

rence  (2.7)  to  approximate  F j (x)  at  the  same  set  of  points, 
using  linear  or  quadratic  interpolation  to  approximate  Ffl(x) 
at  points  x Xj  for  j i a.  Computations  were  performed  with 
a uniform  grid,  i.e.,  x^  * jh,  where  h ■ i/m.  (It  might  be 
more  efficient  to  use  a non-uniform  grid,  because  of  the  log- 
arithmic singularity  of  F^(x)  at  the  origin.)  Using  several 
different  h,  we  found  that  the  error  in  the  computed  value 

of  F (x)  was  0(h),  for  fixed  n and  x.  The  accuracy  could  be 
n 2 

Improved  to  0(h  ) or  better  by  using  Richardson  extrapolation. 
For  example,  using  m ■ 1920  , 3840  and  7680,  we  obtained  Fn(x) 
to  eight  decimal  places  (8D)  for  n £ 20. 

C.  The  Power  Series  Method 

In  Section  3 we  showed  that  F (x)  ■ a (x)lg(x)  + 0 (x) , 

n n n 

where  the  coefficients  a . and  0 in  the  power  series 


a (x)  • ) a ,x^  and  0 (x)  ■ } 0 .x^  satisfy  certain  re- 

n La  n,  j n *->  n,j 

j-0  j-0 

currence  relations.  Thus,  it  is  possible  to  compute  the  co- 
efficients a , and  0 by  working  with  suitably  truncated 
n » J J 

power  series.  To  avoid  numerical  difficulties  it  is  essential 
to  stay  well  within  the  radius  of  convergence  of  each  series, 
which  ensures  that  the  truncated  terms  are  negligible.  This 


|p 


W 
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is  always  possible.  With  Che  series  truncated  after  the 

first  100  terms,  we  computed  Fr(x)  to  12D,  and  the  results 

agreed  with  those  computed  by  the  discretization  method.  The 

value  K ■ 0.705971246102  should  be  correctly  rounded  to  12D. 

Table  7.1  gives  F (x)  to  4D  for  x ■ 0.1 (0.1)0. 9 and 

n»  1(1)5.  It  is  clear  that  the  distributions  F^(x)  converge 

rapidly.  Table  7.2  gives  the  limit  F (x)  to  10D  for  various 

-12 

x.  The  computed  values  of  Ffl(x)  differ  by  less  than  10 
for  all  n 2 20. 

Table  7.3  gives  the  coefficients  a . » 0 . and  § in 

“»J  i J » J 

the  power  series  a M » 3^(x)>  and  5_(x)  “ F 0+x)»  for 
j s 20.  Note  that  the  coefficients  alternate  in  sign,  and 
their  absolute  values  decrease  monotonically,  for  j 2 2. 

The  values  given  in  Tables  7.2  and  7.3  confirm  several 
Identities  which  may  be  derived  theoretically,  for  example: 


2F.<?>  + *. 

<5>  " 2F.<J>  + 2F« 

<£>  + f 

00^3^  + 3» 

H.,,  - -H 

»,2  “ " 

2a  -j » 
»,1 

18B„>2  + 30 

.j  + O0  + 3/ln(2))0oo>1 

- 0. 

Table 

7.1  : Values  of  F 

— n 

(x)  Co 

4D 

X 

F1  (a) 

F2(x) 

F3(x) 

V*) 

f5(x) 

0.1 

0.3329 

0.2871 

0.2772 

0.2753 

0.2750 

0.2 

0.4967 

0.4478 

0.4370 

0.4349 

0,4346 

0.3 

0.6111 

0.5666 

0.5567 

0.5548 

0.5544 

0.4 

0.6989 

0.6611 

0.6526 

0.6510 

0.6507 

0.5 

0.7699 

0.7394 

0.7325 

0.7312 

0.7310 

0.6 

0.8294 

0.8060 

0,8007 

0.7997 

0.7995 

0.7 

0.8805 

0.8637 

0.8599 

0.8592 

0.8590 

0.8 

0.9251 

0.9144 

0.9120 

0.9115 

0.9114 

0.9 

0.9646 

0.9595 

0.9584 

0.9581 

0.9581 

25 


Table  7.2:  Values  of  F (x)  to  10D 


F <x) 
00 


8 

X 

X 

0.2750116116 

1/3 

0.4345648990 

2/3 

0.5544181563 

iA 

0.6507109442 

3/4 

0.7309648721 

1/6 

0.7994844345 

5/6 

0.8590163978 

l/l2 

0.9114387997 

5/12 

0.9580992159 

7/12 

1.0000000000 

ll/l2 

0.5886652481 
0.8400418266 
0.4981238639 
0.8860223000 
0.33708941 90 
0.9275771715 
0.2420627866 
0.6650572783 
0.7887496125 
C. 9653900331 


Table  7.3:  The  Coefficients  n .8  and  F 

“.J  “,j  S“,J 

J 8 f 


For  integers  u and  v,  let  b(u,v)  be  the  number  of  itera- 
tions required  by  the  binary  Euclidean  algorithm  as  described 
in  Section  1 . Let 


(Kv<u<N 
u,v  odd 
and 

(0(N)  - 2B(N)/(LN/2J(LN/2J  - 1)). 

Thus,$(N)  is  Che  average  number  of  iterations  required  for 
distinct,  odd  u and  v less  than  N,  Table  7.4  gives  B(N), 
0(N)  and  A(N)  -0(N)  - (2>( n/2 ) for  N - 23,  24,  215. 

From  the  results  of  Sections  6 and  7,  we  expect  A(N)  to  con- 
verge to  K ■ 0.705971246...  as  N ■*  “.  In  fact,  the  values 
given  in  Table  7.4  satisfy  0 < K - A(N)  < 2 lg(N)/N,  and  give 
the  approximation 

0(N)  « Klg(N)  - 0.93. 

Table  7.4:  Exact  Counts  for  Small  N (algorithm  RS) 


0.6667 


Other  "Binary"  Euclidean  Algorithms 
s the  algorithm  described 
ts  of  the  E 


T 
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For  example,  Harris  [70]  suggested  an  algorithm  which  uses 
both  division  and  right  shifting,  and  requires  less  itera- 
tions than  the  classical  algorithm,  on  the  average.  Yao  and 
Knuth  [75]  considered  the  "subtractive"  Euclidean  algorithm, 
which  requires  neither  shifts  nor  divisions.  In  this  section 
we  analyze  the  "left-shift"  algorithm  (LS)  mentioned  at  the 
end  of  Section  1.  For  positive  integers  u and  v,  even  or 
odd,  the  algorithm  is  as  follows. 

LS  Binary  Algorithm 

LO:  if  u < v then  interchange  u and  v; 

if  u ■ v or  v ■ 0 then  return  u as  the  GCD  and  halt; 
t ♦-  v; 

while  2t  j£  u do  t *-  2 1 ; 

LI : u «-  u - t; 

go  to  LO  . 

The  Interchanging  of  u and  v can  be  avoided  by  duplicat- 
ing some  of  the  code.  The  "while"  loop  merely  shifts  t left 
until  its  leading  one  bit  is  in  the  same  position  as  that  of 
u,  or  one  position  to  the  right  of  it.  This  may  be  done 
with  a floating-point  "normalize"  instruction,  possibly  fol- 
lowed by  one  right  shift. 

Ue  say  that  an  iteration  is  one  execution  of  step  Li . 

The  expected  number  of  iterations  is  given  by  the  following 
theorem. 

Theorem  8.1 

If  integers  u,  v are  chosen  uniformly  and  independently 
in  (0,N],  the  expected  number  of  iterations  of  algorithm  LS 
is  asymptotically  l^lgCN)  as  N ■*  •,  where 

(8.1)  K2  - 12(ln(2)/TT)2  c ~ 0.875837091, 
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(8.2)  c - l PUHgfe^J  , 

J-1 

and  p(j)  is  defined  in  Section  4. 

Proof 

We  shall  only  sketch  the  proof.  Suppose  u>  v > 0 and 
we  perform  one  iteration  of  the  classical  Euclidean  algo- 
rithm, i.e.,  we  find  q ■ |u/vj,  r * u-qv,  8et  u 4-  v and 
v ♦-  r.  Then  the  new  values  of  u and  v would  be  obtained 
after  exactly  p(q)  iterations  of  algorithm  LS.  [Let 


J-1 


where  m^  > m2  > ...  > a 0.  If  1 * j £ p(q)  , then  the 

j-th  execution  of  step  LI  of  algorithm  LS  replaces  the  cur- 

m 

rent  u by  u-t,  where  t ■ 2 Jv.] 

Let  the  regular  continued  fraction  for  u/v  be 

(8.3)  u/v  - qQ  + l/q1  + 1 / ...  + l/qfel 

so  the  classical  algorithm  requires  k+1  iterations.  From 

k 

the  above  discussion,  algorithm  LS  requires  ^ p(q  ) itera- 

u J 
j-0 

tions  (actually  one  less  if  q^  ■ 1,  because  of  our  test 
"if  u - v ..."). 

Let  E2(N)  be  the  expected  number  of  iterations  for 
algorithm  LS,  and  Ec(N)  be  the  expected  number  for  the  class- 
ical algorithm.  Thus, 


(8.4)  lim  E,(N)/E  (N)  - lira  lim  p(qn> , 

N-*®  1 n-*«  N-»® 

where  p(qn>  Is  the  expected  value  of  p(qn> . From  results 
like  those  of  Khinchin  [35a,  35b,  36], 

(8.5)  lim  lim  p(qfl)  - c, 
n-*«  N-*® 

where  c is  given  by  (8.2).  [Intuitively,  the  probability 
that  qn  - q is  about  IsRgJ^J.  from  <2*9)-]  Al8°’ 

(8.6)  Ec(N)  - 12(ln(2)/rr)2lg(N) 

as  N -*  ® (see  Knuth  [69]).  Thus,  the  result  follows  from 
(8.4). 

The  constant  c is  difficult  to  evaluate  numerically 
from  (8.2).  The  following  lemma  is  much  better  for  numerical 
purposes.  Using  (8.8),  we  found 

(8.7)  c = 1.49930818096 
very  easily. 

Lemma  8.1 

If  c is  defined  by  (8.2),  then 

i • 

J j 

• i 

(8.8)  c - 2 + ^ igrti+2--1)  | j 
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(8.10) 


1 + 


1 


21n(2) 


ln(rr) 


“ V + 2 


" (-i?JCOL  . 

^2  J2J(2J-1)J 


Here,  v “ 0.5772...  is  Euler's  constant,  T(x)  is  the  Gamma 
function,  and  £(j)  is  the  Riemann  Zeta  function. 


Sketch  of  Proof 

Splitting  the  sum  in  (8.2)  into  odd  and  even  indices, 
and  using  p(2j+1)  - p(j)  + 1 and  p(2j)  - p(j),  gives 


(8.11)  c 


I 

j-0 


lg 


1+1 7(2 1+1.) 

l+l/(2j+2) 


^ P(j)lg 
j-1 


itl M.  . 

1+1/ (2  j+2) 


Continuing  the  splitting  process  eventually  gives 


(8.12)  c 


00  00 

I Lu 


l+l/  (2k(  j4y)  ) 
l+l/(2k(j+l)) 


From  Stirling's  approximation, 

(8.13)  "ff  [l+x/(j+y)l  ~ nXT(y)/r(x+y) 

J-0 


as  n - so  (8.12)  gives 


(8.14)  c 


lg 


r(j)ni+2"k) 

r(i+2'k) 


From  the  well-known  identity 

(8.15)  r(x)T(x+j)  - T(2x)T(j)21  2x 


(8.17)  lnr(1+x)  - /lnr(n+x)\ 


^ ln(1+x/k) 

k-1 


(8.18) 


lim 

n-*» 


“ n-1 

xln(n)  - ^ ln(1+x/k) 
k-1 


(8.19)  - -yx  + £ (-x)J  . 

j“2 

(8.9)  follows  from  (8.8)  by  putting  x - 2 ^ in  (8.19)  and 
summing  over  k ■ 1,  2,  ...  . The  proof  of  (8.10)  is  similar. 

Numerical  Results  for  Algorithm  LS 

For  integers  u and  v,  let  82(11, v)  be  the  number  of  itera- 
tions required  by  algorithm  LS, 

(8.20)  B2(N)  - ^ b2(u,v), 

0<v<usW 

(8.21)  (3zm  - 2B2(N)/[N(N-1)], 

and 


\ 

\ 


- ..iM  1.1  - .,..1  . 
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(8.22)  Aj  (N)  -^(N)  -<#2( N/2). 

Table  8.1  gives  B2 (N) , ^ (N)  and  ^ (N)  for  N « 22,  23,...,212 
(compare  Table  7.4  for  algorithm  RS). 

Table  8.1 : Exact  Counts  for  Small  N (algorithm  LS) 


N 

B2(N) 

£2(N) 

^(N) 

23 

8 

1.3333 

0.3333 

2? 

55 

1.9643 

0.6310 

2c 

305 

2.5417 

0.5774 

2g 

1625 

3.2762 

0.7345 

2y 

8135 

4.0352 

0.7590 

28 

39282 

4.8329 

0.7977 

184670 

5.6578 

0.8249 

210 

217 

12 

851566 

6.5096 

0.8519 

3860856 

7.3712 

0.8615 

17268497 

8.2383 

0.8671 

76392955 

9.1090 

0.8707 

From  Theorem  8.1,  we  expect 

(8.23)  lim  ^(N)  - K2  « 0.875837, 

If1**0 

and  the  numerical  results  support  this  prediction. 


Summary 

Table  8.2  summarizes  the  average  and  worst-case  behavior 
of  four  algorithms:  the  classical  algorithm,  the  RS  and  LS 

binary  algorithms,  and  the  subtractive  algorithm  of  Yao  and 
Knuth  (75].  The  subtractive  algorithm  is  of  theoretical 
interest  only.  The  choice  of  which  of  the  other  three  algo- 
rithms is  to  be  preferred  depends  on  the  instruction  set  and 
instruction  timing  of  the  machine  used. 
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Table  8.2:  Comparison  of  Various  Euclidean  GCO  Algorithms 


Algorithm 

★ 

Average  iterations 

★ 

Maximum  Iterations 

Classical 

0.5842lg(N) 

1 ,4404 ig(N) 

RS  Binary 

0.70601g(N) 

lg(N) 

LS  Binary 

0.8758lg(N) 

1 .4404 lg(N) 

Subtractive 

0.2921 (lg(N))2 

N 

Notes:  1.  Lover  order  terms  are  neglected  (in  most  cases 

they  are  0(1 ))  . 

2.  An  iteration  of  one  algorithm  (e.g.,  the  binary 
algorithm)  may  take  less  time  than  an  iteration 
of  another  algorithm  (e.g.,  the  classical  algo- 
ri thm) . 
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